A Boussinesq system for two-way propagation 

of interfacial waves 



Hai Yen Nguyen, Frederic Dias a 

a CMLA, ENS Cachan, CNRS, PRES UniverSud, 61, avenue du President Wilson, 

94230 Cachan cedex, France 



Abstract 

The theory of internal waves between two layers of immiscible fluids is important 
both for its applications in oceanography and engineering, and as a source of inter- 
esting mathematical model equations that exhibit nonlinearity and dispersion. A 
Boussinesq system for two-way propagation of interfacial waves in a rigid lid con- 
figuration is derived. In most cases, the nonlinearity is quadratic. However, when 
the square of the depth ratio is close to the density ratio, the coefficients of the 
quadratic nonlinearities become small and cubic nonlinearities must be considered. 
The propagation as well as the collision of solitary waves and/or fronts is studied 
numerically. 



1 Introduction 



As emphasized by Helfrich & Melville [20] in their recent survey article on long 
nonlinear internal waves, observations over the past four decades have demon- 
strated that internal solitary-like waves are ubiquitous features of coastal 
oceans and marginal seas. Solitary waves are long nonlinear waves consist- 
ing of a localized central core and a decaying tail. They arise whenever there 
is a balance between dispersion and nonlinearity. They have been proved to 
exist in specific parameter regimes, and are often conveniently modelled by 
Korteweg-de Vries (KdV) equations or Boussinesq systems. As explained by 
Evans & Ford [16], the differences between "free-surface" and "rigid lid" inter- 
nal waves are small for internal waves of interest. Therefore the "rigid lid" con- 
figuration remains popular for investigating internal waves even if it does not 
allow for generalized solitary waves, which are long nonlinear waves consisting 
of a localized central core and periodic non-decaying oscillations extending to 
infinity. Such waves arise whenever there is a resonance between a linear long 
wave speed of one wave mode in the system and a linear short wave speed of 
another mode [17]. 
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When dealing with interfacial waves with rigid boundaries in the framework 
of the full Euler equations, the amplitude of the central core is bounded by 
the configuration. In the case of solitary waves, it is known that when the 
wave speed approaches a critical value the solution reaches a maximum am- 
plitude while becoming indefinitely wider; these waves are often called 'table- 
top' waves. In the limit as the width of the central core becomes infinite, the 
wave becomes a front [13]. Such behavior is conveniently modelled by an ex- 
tended Korteweg-de Vries (eKdV) equation, i.e. a KdV equation with a cubic 
nonlinear term [18]. Sometimes the terminology 'modified KdV equation' or 
'Gardner equation' is also used. KdV-type equations only describe one-way 
wave propagation. The natural extension toward two-way wave propagation 
is the class of Boussinesq systems. We will derive two sets of Boussinesq sys- 
tems, one with quadratic nonlinearities and another one with quadratic and 
cubic nonlinearities. We will use the terminology 'extended' for a Boussinesq 
system with both quadratic and cubic terms. Some questions arise when deal- 
ing with 'table-top' solitary waves. What are their properties? How do they 
interact? The main goal of this work is to learn more about these waves by 
studying and integrating numerically an extended Boussinesq system which al- 
lows a comparison between fronts and the more standard solitary waves. More 
general models have also been derived by Choi & Camassa [9]. They consid- 
ered shallow water as well as deep water configurations. In the shallow water 
case, their set of equations is the two-layer version of the Green-Naghdi equa- 
tions. The equations derived in [9] were recently extended to the free-surface 
configuration [2]. Solitary waves for two-layer flows have also been computed 
numerically as solutions to the full incompressible Euler equations in the pres- 
ence of an interface by various authors - see for example [22]. Similarly fronts 
have been computed for example in [13,14]. 

The paper is organized as follows. In §2, we present the governing equations 
and the corresponding boundary conditions. A first Boussinesq system of three 
equations is derived in §3. Then it is shown in §4 how to reduce this system 
to a system of two equations, one for the evolution of the interface shape and 
the other one for the evolution of a combination of the horizontal velocities in 
each layer. The numerical scheme and the numerical solutions are described 
in §5. Results are shown for the propagation of a single wave, for the co- 
propagation of two waves and for the collision of two waves of equal as well 
as unequal sizes. When the square of the depth ratio is close to the density 
ratio, the coefficients of the quadratic nonlinearities become small and cubic 
nonlinearities must be considered. An extended Boussinesq system is derived 
in §6. Numerical solutions of the extended Boussinesq system are described 
in § 7. In particular, the collision of 'table-top' waves is considered. A short 
conclusion is given in § 8. In the Appendices, we provide very accurate results 
for wave run-up and phase shift, as well as some intermediate steps in the 
derivation of the extended Boussinesq system. 
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(a) in physical space 



(b) in dimensionless variables 
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Fig. 1. Sketch of solitary waves propagating at the interface between two fluid layers 
with different densities p' and p. The top and the bottom of the fluid domain are 
flat and rigid boundaries, located respectively at z* = h! and z* = —h. (a) Sketch 
of a solitary wave of depression in physical space; (b) Sketch of a solitary wave of 
elevation in dimensionless coordinates, with the thickness h of the bottom layer 
taken as unit length and the long wave speed c as unit velocity. The dashed lines 
represent arbitrary fluid levels 9 and 1 + H — 9' in each layer. The dimensionless 
number H is equal to h'/h. 

2 Governing equations 

The origin of the systems of partial differential equations that will be derived 
below is explained in this section. The methods are standard, but to our 
knowledge some of these equations are derived for the first time. 

Waves at the interface between two fluids are considered. The bottom as well 
as the upper boundary are assumed to be flat and rigid. A sketch is given in 
Figure 1. The analysis is restricted to two-dimensional flows. In other words, 
there is only one horizontal direction, x*, in addition to the vertical direction, 
z* . The interface is described by z* = r]*(x*,t*). The bottom layer Q t * = 
{(x*,z*) : x* G K, — h < z* < i]*(x*,t*)} and the upper layer Q! t * = {(x*,z*) : 
x* G M.,r]*(x* ,t*) < z* < h'} are filled with inviscid, incompressible fluids, 
with densities p and p' respectively. All quantities related to the upper layer 
are denoted with a prime. All physical variables are denoted with a star. 

In addition the flows are assumed to be irrotational. Therefore we are dealing 
with potential flows and only stable configurations with p > p' are considered. 
Velocity potentials <p* — </>* ((x* , z*) , t*) in Q t * and <jf' = <f>*'((x*, z*), t*) in Q' t * 
are introduced, so that the velocity vectors v* and v*' are given by 



v* = V0*, 
v*' = V0*'. 




(2) 
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Writing the continuity equations in each layer leads to 



</W+ 0k* = for -h<z* <r)*(x*,t*), (3) 
0*^* + 0:'*,* = for rj*(x*, t*) < z* < ti. (4) 

The boundary of the system {£l t *,£l' t *} has two parts: the flat bottom z* = —h 
and the flat roof z* = ti. The impermeability conditions along these rigid 
boundaries give 



0** = at z* = -h, (5) 

0*1 = at / = ti. (6) 

The kinematic conditions along the interface, namely D(i]* — z*)/Dt* = 0, 
give 



77* = 0^ - 0*77* at z* = r]*(x*,t*), (7) 
77; = 0:1 - 0;'^ at z* = r]*(x*,t*)- (8) 

The dynamic boundary condition imposed on the interface, namely the con- 
tinuity of pressure since surface tension effects are neglected, gives 



where g is the acceleration due to gravity. The system of seven equations (3)- 
(9) represents the starting model for the study of wave propagation at the 
interface between two fluids. Combined with initial conditions or periodicity 
conditions, it is the classical interfacial wave problem, which has been studied 
for more than a century. A nice feature of this formulation is that the pressures 
in both layers have been removed. In some cases, it is advantageous to keep 
the pressures in the equations. For example, Bridges & Donaldson [8] in their 
study of the criticality of two-layer flows provide an appendix on the inclusion 
of the lid pressure in the calculation of uniform flows. In the next sections, we 
will derive simplified models based on certain additional assumptions on wave 
amplitude, wavelength and fluid depth. 



3 System of three equations in the limit of long, weakly dispersive 
waves 



The derivation follows closely that of [5] for a single layer. Let us now consider 
waves whose typical amplitude, A, is small compared to the depth of the 
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bottom layer h, and whose typical wavelength, £, is large compared to the 
depth of the bottom layeiGl Let us define the three following dimensionless 
numbers, with their characteristic magnitude: 




a At 2 



1. 



Here S is the Stokes number. Let us also introduce the dimensionless density 
ratio r as well as the depth ratio H: 



Obviously r takes values between and 1, the case r = corresponding to 
water waveg[f] while the case r »s 1 corresponds to two fluids with almost the 
same density such as an upper, warmer layer extending down to the interface 
with a colder, more saline layer. The depth ratio takes theoretical values be- 
tween and oo but as said above values H <C 1 or H 3> 1 should be avoided 
in the framework of our weakly nonlinear analysis. 

The procedure is most transparent when working with the variables scaled in 
such a way that the dependent quantities appearing in the problem are all of 
order one, while the assumptions about small amplitude and long wavelength 
appear explicitly connected with small parameters in the equations of motion. 
Such consideration leads to the scaled, dimensionless variables 

x * = £ x , z * = h(z-l), rf = Ar), t* = £t/c , <p* = gAi^/co, 0*' = gA£<j>'/cQ, 

where Co = \fgh- The speed c , which represents the long wave speed in the 
limit r — > 0, is not necessarily the most natural choice for interfacial waves. 
The natural choice would be to take 



which is the speed of long waves in the configuration shown in Figure 1. It 
does not matter for the asymptotic expansions to be performed later. 



There is some arbitrariness in this choice since there are two fluid depths in the 
problem. We could have also chosen the depth of the top layer as reference depth. 
In fact, we implicitly make the assumption that the ratio of liquid depths is neither 
too small nor too large, without going into mathematical details. Models valid for 
arbitrary depth ratio have been derived for example by Choi & Camassa [9]. 
2 In a recent paper, Kataoka [21] showed that when H is near unity, the stability of 
solitary waves changes drastically for small density ratios r. Therefore one must be 
careful in evaluating the stability of air-water solitary waves. In other words, there 
may be differences between r = and the true value r = 0.0013. 
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In these new variables, the set of equations (3)-(9) becomes after reordering 



+ (Pzz 


= 


in 


< z < 1 + ar], 


(10) 




= 


on 


z = 0, 


(11) 




= 


on 


z = 1 + ar], 


(12) 




= 


in 


l + ar]<z<l + H, 


(13) 




= 


on 


z=l + H, 


(14) 




= 


on 


z = 1 + ar], 


(15) 



+ <t>t + + = r (v + <t>' t + + l^'^j on z = 1 + 077.(16) 



We represent the potential as a formal expansion, 

00 

<K{x,z),t)= , £f m (x,t)z m . 



m=0 



Demanding that <fi formally satisfy Laplace's equation (10) leads to the recur- 
rence relation 

(m + 2)(m+l)f m+2 (x,t) = -/3(f m (x,t)) xx , Vm = 0, 1, 2, . . . . (17) 



Let F(x,t) = fo(x,t) denote the velocity potential at the bottom z = and 
use (17) repeatedly to obtain 

Equation (11) implies that fi(x,t) = 0, so 

f 2 k+i(x,t) = 0, Vk = 0,1,2,..., (18) 



and therefore 



, (1 (2k)\ dx 2k 

Let dF(x,t)/dx = u(x,t). Substitute the latter representation into (12) to 
obtain 

r] t + u x + a(ur]) x - ^/3u xxx - ^a(3(r)u xx ) x + ^(3 2 u xxxxx + 0(/3 3 ) = 0.(19) 
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Similarly we represent the potential <p' as a formal expansion, 

CO 

<tf((x,z),t)= "£f^X,t)(l + H-z) m . 



m=0 



Demanding that 0' formally satisfy Laplace's equation (13) leads to the re- 
currence relation 

(m + 2)(m + l)fi +2 (x, t) = -P{f' m {x, t)) xx , Vm = 0, 1, 2, . . . . (20) 



Let F'(x, t) = f' Q (x, t) denote the velocity potential on the roof z — 1 + H and 
use (20) repeatedly to obtain 

nit) _t^?!M VA;-0 12 

f (xt) (-W^ffat) V A;-0 12 

Equation (14) implies that f[(x,t) = 0, so 

& k+1 (x,t) = 0, Vk = 0,1,2,..., (21) 



and therefore 



\2k 



k-- 



" (2k)\ dx 2k 



Let dF'{x,t)/dx = u'{x,t). Substitute the latter representation into (15) to 
obtain 

1 1 

rj t - Hu' x + a(u'rj) x + -f3H 3 u' xxx - -a(3H 2 (r]u' xx ) x 

-^(3 2 H 5 u xxxxx + O((3 3 )=0. (22) 
It is important at this stage that H = 0(1). 

Substitute the representations for <fi and 0' into the dynamic condition (16) to 
obtain the third equation 



(1 _ r)r] + Ft _ rF ' t _ \p (u xt _ rH \'^ 
-a(3r](u xt + rHu' xt ) + ^-(3 2 (u xxxt - rH 4 u' xxxt ) 
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+ ^a{u 2 - I3uu xx ) - ^ar{u' 2 - (3H 2 u'u' xx ) + ^a(3{u 2 x - rH 2 u' 2 ) + 0((3 3 ) = 0. 
Differentiating with respect to x yields 

(f - r)r) x + ut- ru' t - ^[3{u xxt - rH 2 u' xxt ) + a{uu x - ru'u' x ) 
-ap(r]u xt ) x - a(3rH(r]u' xt ) x + —(3 2 {u xxxxt - rH 4 u' xxxxt ) 
- l -afi{uu xx - rH 2 u'u' xx ) x + a(3{u x u xx - rH 2 u' x u' xx ) + 0{(3 3 ) = 0. (23) 

The three equations (19), (22) and (23) provide a Boussinesq system of equa- 
tions describing waves at the interface r](x, t) between two fluid layers based on 
the horizontal velocities u and u' along the bottom and the roof, respectively. 
It is correct up to second order in a, (5. 

One can derive a class of systems which are formally equivalent to the system 
we just derived. This will be accomplished by considering changes in the de- 
pendent variables and by making use of lower-order relations in higher-order 
terms. Toward this goal, begin by letting w(x,t) be the scaled horizontal ve- 
locity corresponding to the physical depth (1 — 9)h below the unperturbed 
interface, and w'(x,t) be the scaled horizontal velocity corresponding to the 
physical depth (H — 6')h above the unperturbed interface. The ranges for the 
parameters and & are < 9 < 1 and < & < H. Note that (9, 9') = (0, 0) 
leads to w = u and w' = u', while (9, 9') = (1, H) leads to both velocities eval- 
uated along the interface. A formal use of Taylor's formula with remainder 
shows that 

w = <f> x \ z= g = (F x - l -?F xxx 9 2 + ^P 2 9 4 F xxxxx ^j + 0(P 3 ) 
= u- \i59 2 u xx + ^[3 2 9 4 u xxxx + 0(/3 3 ) 
as (3 — > 0. In Fourier space, the latter relationship may be written as 

w=(l + -[39 2 k 2 + —B 2 9 4 k A \ u + 0(p 3 ). 
\ 2 24 / 

Inverting the positive Fourier multiplier yields 



1 



1 



u - 



l + + w + 0(JP) 

1 _ l {3e 2 k 2 + -6 2 9 4 k 4 ) w + 0(P 3 ) 
v 2 24 / 
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as (3 — > 0. Thus there appears the relationship 

u = w + ^(36 2 w xx + ^(3 2 6 4 w xxxx + 0{f3 3 ). (24) 



Similarly 

w' = tiU+H-r = (K - + ^ 2 F> XXXX 0*) + 0((3 3 ) 

= «' - \(3d' 2 u' xx + ±(3 2 9' 4 u' xxxx + 0{(3 3 ) 

and 

w> = (i + l -pe' 2 k 2 + Y A P 2 e"k^ u> + o((3 3 ). 

Inverting the positive Fourier multiplier yields 

u'=(l- -f36' 2 k 2 + —pe'^) w' + 0{(3 3 ) 
V 2 24 / 

and thus the relationship 

u' = w> + \pe' 2 w' xx + ^ 2 ,A w' xxxx + 0((3 3 ). (25) 



Substitute the expressions (24) and (25) for u and u' into (19) and (22), 
respectively, to obtain 



Vt + Wx + a(wr]) x + ^(3 (o 2 - ^ 



w 7 



+i a /3(0 2 - l)( V w xx ) x + ^-B 2 (e 2 - l -]w xxxxx + O((3 3 ) = 



2 / v /v i ■ 24 V 5 



jfe - + a(«; , 77) a - (V 2 - \h 2 \ ,r\ 



(26) 



+\a(3 (6> 2 - H 2 ) {r)w' xx ) x - ^(3 2 H (> - \h 2 ^ w' xxxxx + O((3 3 ) = 0. 

(27) 

Substitute the expressions (24) and (25) for u and u' into (23) to obtain 



(1 — r)i] x + w t — rw' t H — (3 (6 2 — l)w — r(0' 2 — H 2 )w' + a(ww x — rw'w' x ) 



2" 



xxt 



1 ,2 
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[{r]w xt ) x + rH(r]w' xt ) x ] + -af3 [(9 2 - l)ww xxx - r(6' 2 - H 2 )w'w' xxx 

+ l -aP [(9 2 + l)w x w xx - r(9' 2 + H 2 )w' x w' xx ] + 0((3 3 ) = 0. 

(28) 



The system of three equations (26)-(28) is formally equivalent to the previous 
system but it allows one to choose the fluid levels 9 and 9' as reference for the 
horizontal velocities. Among all these systems that model the same physical 
problem one can select those with the best dispersion relations. Neglecting 
terms of 0(a 2 , (3 2 , a(3), the system (26)-(28) reduces to 

rj t + w x + a(wrj) x + \i5{9 2 - \)w xxx = 
Vt - Hw' x + a(w' V ) x - \PH{9> 2 - \H 2 )w' xxx = (29) 
(1 — r)i] x + w t — rw' t + \l3[{9 2 — l)w — r(9' 2 — H 2 )w'] xxt + a(ww x — rw'w' x ) = 

4 System of two equations 

The systems obtained in the previous section are not appropriate for numerical 
computations. One would like to obtain a system of two evolution equations 
for the variables r] and W = w — rw'. In fact, Benjamin and Bridges [3] (see 
also [12,11,1] ) formulated the interfacial wave problem using Hamiltonian 
formalism and showed that the canonical variables for interfacial waves are 



rj*(x*, t*) and p0*(x*, rj*, t*) - p'(p* (x*, i]*, t*). 



At leading order, the first two equations of system (29) give 




Assuming the fluids to be at rest as x — > oo, one has w = —Hw'. Therefore 



H 



W + 0((3), w' = 



-1 



w + o(p). 



(30) 



w = 



r + H 



r + H 



Adding H times the first equation to r times the second equation of system 
(29) yields 



(r + H)i] t + H{w — rw') x + a[(Hw + rw')rf\ 



(31) 



+f p{(9 2 - \)w xxx - t(9' 2 - \H 2 )w' xxx ] = 0. 
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Using (30) and neglecting higher-order terms, one obtains 
where 

S = (9 2 -l) + ^-(e' 2 -H 2 ). 

In the third equation of system (29), the term with the xxt— derivatives can 
be written as 



1 

2' 



/? (tf 2 - I)*;** - r{6' 2 - H 2 )w' xxt ] = l(3^W xxt 



The quadratic terms of the third equation of system (29) can be written as 



H 2 — r 

a(ww x - rw'w' x ) = a h)2 WWx ' 



Then the third equation of system (29) becomes 

W t = -(1 - r)rj x - \p^W xxt - a^y 2 WW x . 

The final system of two equations for interfacial waves in the limit of long, 
weakly dispersive waves, can be written in terms of the horizontal velocities 
at arbitrary fluid levels as (in dimensionless form) 



(32) 



(33) 



Vt = -4hW x - af^(W V ) x - P (1^ + iZftgp) W X3 
( W t = -(1 - r) Vx - a^=fyWW x - \(5^W xxU 

or as (in physical variables) 

= -hd x W*, - d A (W*r]*) x * - h 3 d 2 W** x * x *, 

w * = _ g (! _ r)r} .. _ d 4 w*w** - h 2 d 3 w; w , 

where 

Notice that Choi & Camassa [9] also derived a system of two equations (see 
their equations (3.33) and (3.34)), but it is different from ours. In particular, 
their coefficient d 2 is equal to 0, and their equation for W t possesses an extra 
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quadratic term r\r\ x . The reason is that their ''W 1 is the mean horizontal velocity 
through the upper layer. The value of S which best approximates the Choi 
& Camassa equations is S = — 1(1 + rH). Indeed the coefficient d 2 then 
vanishes. This particular value for S can be explained as follows. The leading 
order correction to the horizontal velocity is given by 

w(z) = u - \(3z 2 u xx . 

The value of z, say z = 9, for which the mean velocity 



i 

T j w( dr 




is equal to w(9) is given by 9 = l/y/3. Similarly, one finds 9' = {l/y/S)H for 
the upper layer. Therefore S — — 1(1 + rH). 

Recall that the scaling that led to our Boussinesq system is given by 

X* X t* _ t T]* _ W* _ ^ 

h v//7' h/co yf]3' h aT] ' gh/c 

with Co = \fgh, a < 1, /3 < 1 and a = 0(/3). Linearizing system (33) and 
looking for solutions (r]*,W*) proportional to exp(ikx* — iuot*) leads to the 
dispersion relation 

_ gh{l-r){d l -d 2 k 2 h 2 ) 
F ~ 1 - d 3 k 2 h 2 ' 

Plots of the dispersion relation are given in the next section. Since < 9 < 1 
and < 9' < H, the definition of S implies that 

-1 -rH < S < 0. 

It follows that g?3 < and therefore the denominator 1 — d 3 h 2 k 2 is positive. 
In order to have well-posedness (that is uj 2 /k 2 positive for all values of k), d 2 
must be negative, which is the case if S < — 1(1 + rH). Finally the condition 
we want to impose on S is that 

-(l + rH)<S <--(l + rH). (35) 



It is satisfied if one takes the horizontal velocities on the bottom and on the 
roof (S = —(1 + rH)) or the mean horizontal velocities in the bottom and 
upper layers (S = — 1(1 + rH)), but it is not if one takes the horizontal 
velocities along the interface (S — 0). 
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5 The numerical scheme and numerical solutions 



In order to integrate numerically the Boussinesq system (33), we introduce a 
slightly different change of variables, where the stars still denote the physical 
variables and no new notation is introduced for the dimensionless variables: 

x * V* c „ Tir W* . , 2 H(l-r) 
x = — , V = -T, t= T t*, W = , with c 2 = gh- 



h h h c r + H 

The system (33) becomes 



r] t = -diW x - d 4 {Wr]) x - d 2 W xxx 

1 , (36) 

W t = -~rVx - d A WW x - d 3 W xxt 

di 



with dispersion relation 
u 2 d\ — d 2 k 2 

¥ = di(i-d 3 fc 2 )" 



(37) 



As k — > 0, — > 1. As A; — > oo, 

^ ^ 2(l + rff) 

A; 2 ^ didg + 35 ' 

Typical plots of the dispersion relation (37) are given in Figure 2. Comparisons 
between the approximate and the exact dispersion relations, given by 

oj 2 tanh k tanh kH 

k 2 (^(tanh kH + rtanh k) 

are also shown. A very good agreement is found for small k. Taking the Fourier 
transform of the system (36) gives 

rj t = (d 2 k 2 - d x )ikW - d 4 ikWr] 

= - dl (i-d 3 k 2 ) ikf] - w^dwf ™ 2 ' 

The system of differential equations is solved by a pseudo-spectral method in 
space with a number iV of Fourier modes on a periodic domain of length L. For 
most applications, iV = 1024 was found to be sufficient. The time integration 
is performed using the classical fourth-order explicit Runge-Kutta scheme. 
The time step At was optimized through a trial and error process and was 
found to have a dependence in 1/N. 
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(a) (b) 




0.5 1 1.5 2 0.5 1 1.5 2 



k k 

Fig. 2. Dispersion relation (37) for the Boussinesq system (36) with S = —1 — rH, 
r = 0.9: (a) H = 1.2, (b) H = 0.8. The dashed curves represent the dispersion rela- 
tion for the linearized interfacial wave equations, without the long wave assumption 
(see for example [22]). 



Since the main goal is to study the propagation and the collision of solitary 
waves, we first look for solitary wave solutions of the system (36). As op- 
posed to the KdV equation, there are no explicit solitary wave solutions of the 
Boussinesq system that are physically relevant. Therefore we look for an ap- 
proximate solitary wave solution to (36) as in [4] (see also [15] for the existence 
of solitary wave solutions). The leading-order terms give 

di 

A solution representing a right-running wave is 

W(x-t) = ^-r](x-t). 

Let us look for solutions of system (36) in the form 

W(x,t) = -[ V (x,t) + M(x,t)}, 

where M is assumed to be small compared to r\ and W. Substituting the 
expression for W into (36) and neglecting higher-order terms yields 

Vt = ~Vx ~ M x - ^(r] 2 ) x - 

h _ d ^ . (38) 



Vt = -Vx -M t - ^j-(v 2 )x - d 3 7] xxt 
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Assuming that the solitary wave goes to the right, one has M t —M x . There- 
fore 

Af x = ---r{V )x ~ ^-rVxxx + 7i d 3Vxxt- 
4 d\ zdi 2 

Substituting the expression for M x into one of the equations of system (38) 
yields 

3g?4 , o\ d 2 d$ /nn\ 

Vt + r]x + 4d~i ' x + 2d? xxx + ~2 Vxxt = ^ ' 



This is essentially the model equation that was studied in [7]. 
Looking for solitary wave solutions of (39) in the form 

i] = i] sech 2 [/t(:r + x - Vt)} (40) 

leads to two equations for k and V: 

-V + 1 + 2(d 2 /d l )K 2 - 2d 3 K 2 V = 
d 4 r) - 4:d 2 n 2 + Adid?,K 2 V = 

Solving for k 1 and V yields 

2 d 4 r] , r , . d 4 r] 

K = — 7 V' V = I + , 

4 (d 2 - ditfe - ^d 3 d 4 r)o) 2d i 

and, assuming M(±oo) = 0, one obtains explicitly the following expression 
for M: 

M = ~ 2d?™ + y^- 

For a given pair (r, if), one must only consider values of r]Q which are such that 
k 2 > 0. In addition one has the condition (35) on S. The sign of d 4 depends 
on the relation between H 2 and r. Let us assume first that H 2 > r so that 
d 4 > 0. In order for the condition k 2 > to be satisfied, one needs 



Vo (d>2 - did 3 - ^d 3 d 4 i]Q^ > 0. 



The values of rj for which the left-hand side of the inequality vanishes are 

_ AH(r + H)(1 + rH) 
^701-0, Vo2 - 3{H2 _ ^ . 
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Since S < 0, ^02 < and therefore 7702 < Voi- The coefficient of i] 2 in the 
inequality is positive. Consequently one must have 

AH(r + H)(l + rH) 
Vo > V01 = or 770 < V02 = s(H 2 -r)S ' 

This second branch is not acceptable since 

4H(r + H)(l + rH) 

3(/P-r) > >^rH>-S>0. 



Therefore 



AH(r + H)(l + rH) 

- < -1, 



3(H 2 -r)S 

which gives an amplitude larger than the depth! 

Similarly, when H 2 < r one finds a second branch which is not acceptable. 
The summary of acceptable values for tjq is given in the table 



H 2 - r > 


< 770 < H 


H 2 - r < 


-1 < T] < 



For a "thick" upper layer (H 2 > r), the solitary waves are of elevation, while 
they are of depression for a "thick" bottom layer (H 2 < r). The weakly non- 
linear theory developed in the present section does not provide any bounds 
on the amplitude of the solitary waves. We have added a physical constraint 
based on the fact that both layers are bounded by flat solid boundaries. It 
is well-known in the framework of the full interfacial wave equations (see for 
example [22]) that the rigid top and bottom provide natural bounds on the 
solitary wave amplitudes. As the speed increases, the wave amplitude reaches 
a limit. In the next section, we extend our weakly nonlinear analysis to cubic 
terms so that this effect can be incorporated. 

Once the approximate solitary wave (40) has been obtained, it is possible to 
make it cleaner by iterative filtering. This technique has been used by several 
authors, including [4,6], and is explained in Appendix A. In order to study 
run-ups and phase shifts during collision of solitary waves, it is important to 
use clean solitary waves for the initial conditions. On the other hand, in order 
to show only the qualitative behavior, it is not necessary. Therefore results in 
this Section are given for non-filtered solitary waves. Some results with filtered 
waves are described in Appendix A. 

In Figure 3, we show the propagation of an almost perfect right-running soli- 
tary wave of elevation. Even though all computations are performed with 
dimensionless variables, it is interesting to provide numerical applications for 
a configuration that could be realized in the laboratory [23]. Keeping r = 0.9 
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Fig. 3. An approximate solitary wave propagating to the right. This is a solution 
to the system of quadratic Boussinesq equations (36), with parameters H = 1.1, 
r = 0.9, L = 512, N = 1024, S = -1 - rH, r/ = 0.1. 

as in the figure, one could take for example h — 10 cm, h' — 11 cm (H — 1.1). 
The solitary wave amplitude is 1 cm, its speed c 23.2 cm/s, the length of 
the domain 51.2 m (a bit long!). The plots (b)-(e) would then correspond to 
snapshots at t — 21.5 s, t — 68.9 s, t — 94.8 s and t = 163.7 s. 

In Figure 4, we show the head-on collision of two almost perfect solitary waves 
of elevation of equal amplitude moving in opposite directions. As in the one- 
layer case, the solution rises to an amplitude slightly larger than the sum of 
the amplitudes of the two incident solitary waves (see Appendix A). After the 
collision, two similar waves emerge and return to the form of two separated 
solitary waves. As a result of this collision, the amplitudes of the two result- 
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(f) evolution in time 




Fig. 4. Head-on collision of two approximate solitary waves of elevation of equal 
size. This is a solution to the system of quadratic Boussinesq equations (36), with 
parameters H = 1.2, r = 0.8, L = 512, N = 1024, S = -1 - rH, rf = rf = 0.1, 
where the superscripts £ and r stand for left and right respectively. 



ing solitary waves are slightly smaller than the incident amplitudes and their 
centers are slightly retarded from the trajectories of the incoming centers (see 
again Appendix A). 

In Figure 5, we show the collision of two almost perfect solitary waves of 
depression of unequal amplitudes moving in opposite directions. The numerical 
simulations exhibit a number of the same features that have been observed in 
the symmetric case. 
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(a) t = 
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(f) evolution in time 




Fig. 5. Head-on collision of two almost perfect solitary waves of depression of differ- 
ent sizes. This is a solution to the system of quadratic Boussinesq equations (36), 
with parameters H = 0.6, r = 0.85, L = 256, N = 1024, S = -1 - rH, rf = -0.04, 
t/q = —0.11, where the superscripts £ and r stand for left and right respectively. In 
plot (f), note that —rj(x,t) has been plotted for the sake of clarity. 



In Figure 6, we show the co-propagation of two solitary waves of elevation of 
different amplitudes. A sequence of spatial profiles is shown. The larger one, 
which is faster, eventually passes the smaller one, which is slower. Again there 
is a phase shift after the interaction. The amplitude of the solution r](x, t) never 
exceeds that of the larger solitary wave, nor does it dip below the amplitude 
of the smaller. 
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Fig. 6. Co-propagation of two almost perfect solitary waves of elevation of different 
sizes. This is a solution to the system of quadratic Boussinesq equations (36), with 
parameters H = 1.6, r = 0.95, L = 2 14 , N = 2 14 , S = -1 - rH, r& = 0.1, rf Q = 0.03, 
where the superscripts t and r stand for left and right respectively. 



6 Extended Boussinesq system of two equations with cubic terms 



When \H 2 — r\ is small, one needs to go one step beyond and take into con- 
sideration the cubic terms. Again one would like to obtain a system of two 
equations for the variables rj and W = w — rw'. We derive first a general sys- 
tem of two equations with cubic terms. Then we introduce a specific scaling 
for the case where \H 2 — r\ is small. A lot of terms in the system drop out 
because they are of higher order. 
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The leading order terms lead to the same equation as before, namely w = 
—Hw'. And again 

w = -^—W + 0((3), w ' = ^L-W + 0(P). (41) 
r + H r + H 



At the next order, the first two equations of (29) give 

w x +a{w V ) x +^(3 (o 2 - w xxx = -Hw' x +a{w' V ) x -^(3H (V 2 - ^H 2 ^j w' xxx . 
Since the speeds w and w' vanish as x — > oo one has 

w = _ Hv / + a{w > _ ^ _ \p (h (> _ ^ w ^ + ^ 2 _ ^ 

Using (41) for the terms containing a or/3 and neglecting terms of 0(f3 2 ), one 
obtains 



w = 



w = 



-Hw' 

w 
H 



a 



i + H i (o 12 - Ih 2 ) - (e 2 - f 

r + H 2 r + H 



a 1 + H Wrj+'-P^ 
H(r + H) 1 T 



12 



(e 2 



r + H 



W x 



(42) 
(43) 



In Appendix B, after several substitutions, one obtains the system of two 
equations (B.8) and (B.15). Switching back to the physical variables 



x* = £x, rf = Ar], t* = it I c , W* = gAW/c , with c = y/gh, 
the system (B.8)-(B.15) becomes 

(r + h)t& + hHw;, + ^(w* v %* + itfH ^-^f-^ w;.^. 

+\# (rH(l + H rXnt^ + ^T^ 1 ) W^')- W 



1 l5 f r^((^-l^)-(^-l)) 5 l)2 +rH(g /2_ l g 2)2 \ 

4 U Cr--4-ff12 fi r +ff VV X*X*X*X*X* ~ U > 
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, , 4 ^ rH((e 2 -l)-(e' 2 -H 2 ))((g' 2 -ig 2 )-(e 2 -l)) g (0^l)( 5 e 2 -l)+r(0' 2 -ff 2 )(5e' 2 -g 2 ) 
+n y 4(r+H) 2 2(r+H) 

W* — 



The specific scaling for small values of \H 2 — r|, 

X* X t* t T]* W* - 

T" = 777 — = ~> "77 = ar 7> ~T7 = ^ - r = aC, 

h p h/c a n gh/co 



with c = v / #^> a -C 1, /3 -C 1, a = 0(/5), will lead to a new Boussinesq 
system with cubic terms. A lot of terms in (44)- (45) drop out because they 
are of higher order. Keeping terms of order a 2 and a 4 and going back to 
physical variables, the system of two equations becomes 

'It*— a r+H VV x* 11 \2(r+H) 2 ^3 (r+H) 2 ) VV x*x*x* 

W* = -g(l - r)i£. - \h 2 ^W: w - ^W*W* X * + \ r -^{W^% 



This is the same system as (33) with two extra terms, the cubic terms. We will 
call it a system of extended Boussinesq equations (see also [11]). Linearizing 
(46) gives the same dispersion relation as before. 



7 Numerical solutions of the extended Boussinesq system 



In order to integrate numerically the extended Boussinesq system (46), we 
introduce a slightly different change of variables, where the stars still denote 
the physical variables and no new notation is introduced for the dimensionless 
variables: 

x* rj* c W* 2 ghh'(p-p') ghH(l-r) 

x = — , 7] = — , t — —t , W — , with - 



h' 1 h' h ' c ' p'h + ph' r + H 
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Using the same coefficients as in (34), we rewrite system (46) with the new 
variables as 



rj t = -d x W x - d 2 W xxx - d 4 (Wr]) x + d 5 (Wr] 2 ) x 

(47) 

w t = -(iM)% - d 3 w xxt - d A ww x + d 5 (w 2 v ) x 



where the new coefficient d 5 is equal to 

r(l + H) 2 



d* 



(r + Hf 



When (6, 9') = (0, 0), one recovers the system with horizontal velocities on the 
bottom and on the roof. 

Taking the Fourier transform of the system (47) gives 



r},. = (d 2 k 2 - d 1 )ikW - d 4 ik(Wr]) + d 5 ik(Wr] 2 ), 

(1 - d 3 k 2 )W t = —ikfj - ^ikCW 2 ) + d^ikiW 2 ^). 
d\ 2 

The system of differential equations is integrated numerically with the same 
method as in § 5. 

Again we look for approximate solitary wave solutions to (47). As before we 
look for solutions of the form 

W(x,t) = -[ V (x,t)+M(x,t)}, 

where M is assumed to be small compared to rj and W. Substituting the 
expression for W into (47) and neglecting higher-order terms yields 

1 d& , o\ 1 d 2 1 , , 

M x = -- — [rj ) x - -—r) xxx + -d 3 rj xxt . (48) 



Substituting the expression for M x into one of the equations of system (47) 
yields 

Vt + Vx + j^Vn )x - -jW )x + T^fVxxx + -^Vxxt = 0. (49) 



We have checked that the extended KdV equation (49) is in agreement with 
previously derived eKdV equations such as in [14]. 
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(b) 



Fig. 7. 'Table-top' solitary waves which are approximate solutions of the extended 
Boussinesq system (47). The horizontal velocities are taken on the top and the 
bottom so that S = —(1 + rH). (a) H = 1.8, r = 0.8. The wave speeds V are, going 
from the smallest to the widest solitary wave, V max — V ~ 10~ 3 , 10 -9 , 10~ 15 ; (b) 
H = 0.4, r = 0.9. The wave speeds V are, going from the smallest to the widest 
solitary wave, V max - V ~ 10~ 3 , 10~ 9 , 10" 14 . 

Let V = 1 + Ci be the wave speed, with c\ small. In the moving frame of 
reference X — x — (1 + ci)t, T — t, equation (49) becomes 

-Cir]x+VT+j^(v 2 )x- ( ^-(r] 3 )x+^-Vxxx+ ( ^- [-(1 + c^rjxxx + Vxxt] = 0. 
Looking for stationary solutions and integrating with respect to X yields 

-^+^-f/ + \(fr d >-^h^°- (50) 



Letting 

3 H 2 -r _ r(l + Hf ! H(rH+l) 1 HS 

ai ~2H(r + HY ~ H(r + H) 2 ' 1 ~ 6 r + H 47+# Cl ' 



equation (50) becomes 

-ci?7 + ^iV 2 ~ + ^iVxx = 0. 

It has solitary wave solutions 

v{x) = 1 21) — i-' 3 , with t = iMZpi. 

\Pi) l + ecosh( v /fX) |«i| 
In the fixed frame of reference, the profile of the solitary waves is given by 
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When H 2 > r the solitary waves are of elevation. When H 2 < r they are of 
depression. The parameter e can take values ranging from (infinitely wide 
solution) to 1 (solution of infinitesimal amplitude). Assuming M(±oo) = 0, 
one can compute M explicitly by integrating equation (48) with respect to x: 



Typical approximate solitary waves solutions are shown in Figure 7. Notice 
that the condition \H 2 — r\ small is not really satisfied for the selected values 
of H and r. The reason is that otherwise the waves would have been too 
small to be clearly visible. Of course we still have the conditions on S for 
well-posedness: 



The solitary waves are characterized by wave velocities larger than 1 (ci > 0). 
The maximum wave velocity V max is obtained when e — > 0. One finds c\ — > 



Once the approximate solitary wave (51) has been obtained, it is again possible 
to make it cleaner by iterative filtering. Qualitative results for non-filtered 
solitary waves are given in this Section. Some accurate results for run-ups and 
phase shifts with filtered waves are described in Appendix A. 

In Figure 8, we show the head-on collision of two almost perfect 'table-top' 
solitary waves of elevation of equal amplitude moving in opposite directions. 
As in the case with only quadratic nonlinearities, the solution rises to an 
amplitude larger than the sum of the amplitudes of the two incident solitary 
waves. After the collision, two similar waves emerge and return to the form 
of two separated 'table-top' solitary waves. As a result of this collision, the 
amplitudes of the two resulting solitary waves are slightly smaller than the 
incident amplitudes and their centers are slightly retarded from the trajectories 
of the incoming centers. 

In Figure 9, we show the collision of two almost perfect solitary waves of 
depression of equal amplitude moving in opposite directions. The numerical 
simulations exhibit the same features that have been observed in the elevation 
case. 



M = 





af/6(3i, so that 
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Fig. 8. Head-on collision of two approximate 'table-top' elevation solitary waves of 
equal size. This is a solution to the system of cubic Boussinesq equations (47), with 
parameters H = 0.95, r = 0.8, L = 4096, N = 1024, S = -1-rH, V max -V ~ 10~ 17 . 



In Figure 10, we show the collision of an almost perfect 'table-top' solitary 
wave of elevation with a solitary wave of elevation moving in the opposite 
direction. The numerical simulations exhibit a number of the same features 
that have been observed in the symmetric case. The phase lag is asymmetric, 
with the smaller solitary wave being delayed more significantly than the larger. 



Note that in the quadratic as well as in the cubic cases, it is not possible 
to consider the collision between a solitary wave of depression and a solitary 
wave of elevation. Indeed the sign of H 2 — r determines whether the wave is 
of elevation or of depression. 
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(c) t = 700 




(d) t = 1000 






-0.005 
-0.01 

-0.015 
-0.02 

-0.025 



(e) * = 2200 



(f) evolution in time 




Fig. 9. Head-on collision of two approximate 'table-top' depression solitary waves of 
equal size. This is a solution to the system of cubic Boussinesq equations (47), with 
parameters H = 0.9, r = 0.85, L = 4096, N = 1024, S = -1-rH, V max -V ~ 10~ 14 . 
In plot (f), note that —r](x,t) has been plotted for the sake of clarity. 



8 Conclusion 



In this paper, we derived a system of extended Boussinesq equations in order 
to describe weakly nonlinear waves at the interface between two heavy fluids 
in a 'rigid-lid' configuration. To our knowledge we have described for the 
first time the collision between 'table-top' solitary waves. The extension to 
a 'free-surface' configuration and to arbitrary wave amplitude is left to future 
studies. Indeed, since the waves we considered are only weakly nonlinear, we 
do not have to worry about the resulting wave reaching the roof or the bottom. 
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Fig. 10. Head-on collision of a solitary wave of elevation and of a 'table-top' solitary- 
wave of elevation. This is a solution to the system of cubic Boussinesq equations 
(47), with parameters H = 0.95, r = 0.8, L = 2048, N = 1024, S = -1 - rH, 



V r ~ 10 



-11 



However, in a fully nonlinear regime, this could happen. Indeed the maximum 
amplitude A for 'table-top' solitary waves is given by 

A _ H - y/r 
h ~ 1 + v/F ' 

Take the case where H 2 > r. It is easy to see that while A/h is always smaller 
than H, 2 A/h can exceed H, so that the resulting wave will hit the roof. 
Therefore it will be interesting to consider the collision of solitary waves of 
arbitrary amplitudes by using the full Euler equations. On the other hand, for 
'table-top' solitary waves of depression, the resulting wave cannot touch the 



28 



bottom. 



A Additional results on run-ups and phase shifts 

In this appendix, we provide accurate results on run-ups and phase shifts. 
The terminology 'run-up' denotes the fact that during the collision of two 
counterpropagating solitary waves the wave amplitude increases beyond the 
sum of the two single wave amplitudes. Since run-ups and phase shifts are 
always very small, they must be computed with high accuracy. This is why it is 
important to clean the solitary waves obtained by the approximate expressions 
(40) or (51). We proceed as follows. We begin with an approximate solution, 
let it propagate across the domain, truncate the leading pulse, use it as new 
initial value by translating it to the left of the domain, let it propagate again 
and distance itself from the trailing dispersive tail, truncate again, and repeat 
the whole process over and over until a clean, at least to the eye, solitary wave 
is produced. Then we use this new filtered solution as initial guess to study 
the various collisions. 

For solitary wave solutions to the system of equations with quadratic nonlin- 
earities (36), the behavior is the same as the behavior shown for example in 
[10]. In particular we obtain pictures that look very similar to their Figure 2 
for the phase shift resulting from the head-on collision of two solitary waves 
of equal height, to their Figure 4 for the time evolution of the maximum am- 
plitude of the solution (it rises sharply to more than twice the elevation of the 
incident solitary waves, then descends to below this level after crest detach- 
ment, and finally relaxes back to almost its initial level) and to their Figure 12 
for the asymmetric head-on collision of two solitary waves of different heights. 

Since the main contribution of the present paper is the inclusion of cubic 
terms in addition to the quadratic terms, we focus on results for the extended 
Boussinesq system (46). Figure A.l shows the effect of cleaning. In Figure 
A. 2, the collision between two clean 'table-top' solitary waves (the clean- 
ing has been applied 400 times) is shown. Their speed is V = 1.00183358. 
The amplitude before cleaning was ?? max = 0.063476. After iterative clean- 
ing, it reached r) maK = 0.06812113. The run-up during collision is extremely 
small: indeed r/ max = 0.13624323 at collision, which is slightly larger than 
2 x 0.06812113 = 0.13624226. The phase shift is also very small. In Figure 
A. 3, the collision between the clean 'table-top' solitary wave of the Figure A.l 
and a clean solitary wave (the cleaning has been applied 230 times) is shown. 
The maximum amplitude is greater than the sum of the two wave amplitudes. 
The speed of the smaller wave is V = 1.0015. Its amplitude before cleaning 
was i] max = 0.03647847. After iterative cleaning, it reached r] max = 0.03719492. 
The run-up during collision is again extremely small, even if it is larger than 
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Fig. A.l. Flat solitary wave produced by iterative cleaning. This is a solution to the 
system of extended Boussinesq equations (47). (a) Difference in the profile before 
(solid line) and after (dashed line) cleaning, (b) Profile of the approximate solitary 
wave (51) after one propagation across the domain, (c) Profile (b) after cleaning 
and translation to the left of the domain, (d) Profile after several cleanings. Notice 
the change of scale in the vertical axis, (e) Evolution of the maximum amplitude 
f?max as cleaning is repeated over and over. The amplitude reaches an asymptotic 
level. 



in the previous case: indeed r/ max = 0.10556057 at collision, which is slightly 
larger than 0.06812113 + 0.03719492 = 0.10531605. The phase shift is very 
small and the crest trajectory shows an interesting path. The overall conclu- 
sion is that run-ups and phase shifts are smaller for 'table-top' solitary waves 
than for 'classical' solitary waves. 
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Fig. A. 2. A collision between two clean 'table-top' solitary waves of equal height. 
This is a solution to the system of extended Boussinesq equations (47). (a) Initial 
profiles, (b) Time evolution of the amplitude i] mSLX . (c) Crest trajectory. 

B Intermediate steps in the derivation of the extended Boussinesq 
system with cubic terms 



Adding H times equation (26) to r times equation (27) yields 



(r + H)i] t + H{w x — rw' x ) + a[(Hw + rw')r]\ a 



(0 2 - \)w xxx - r(9' 2 - \h 2 )w' xi 



1 



5H o2 
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Q2 

1 2 



+ -a(3 [H(6 2 - l)( V w xx ) x + r{6> 2 - H 2 )( V w' xx ) x 



5 1 "'' 



= 0. 



(B.l) 



Let us replace the variables w and w' in (B.l) by their expressions (42)-(43) 
in terms of W and let 



1 + H 



G = H 



r + H' 

We consider all the terms one by one: 



(6> 2 - \H 2 ) - (6 2 - |) 
r + H 
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Fig. A. 3. A collision between a clean solitary wave and a clean 'table-top' solitary- 
wave. This is a solution to the system of extended Boussinesq equations (47). (a) 
Initial profiles, (b) Time evolution of the amplitude r/ max . (c) Crest trajectory. 



(1) Term H{w x -rw' x ) 

H(w x - rw' x ) = HW X 

(2) Term a[(Hw + rw')r)] x 
From equation (42), we have 
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so that 
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w - rw' = -(r + H)w' - aFWrj + -f3GW xx , 
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r + H r + H ' 2' r + H 

Similarly from equation (43), we obtain 



w = 



r + H r + H 

Combining (B.3) and (B.4) yields 
H 2 
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Therefore 
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(3) Term in w xxx and w' xxx 

Combining (B.3) and (B.4) yields 
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(4) Term «n {f]w xx ) x and {r)w' xx ) 3 
Using (41) yields 
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(5) Term in w xxxxx and w' x , 
Using (41) yields 
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Combining all terms (B.2)-(B.7) yields the first equation of the extended 
Boussinesq system 
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We proceed the same way for equation (28). 
(1) Term in w xx t and w' xxt 
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Combining all terms (B.9)-(B.14) yields 
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